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We generalize the hydrodynamic lattice gas model to include arbitrary numbers of par- 
ticles moving in each lattice direction. For this generalization we derive the equilibrium 
distribution function and the hydrodynamic equations, including the equation of state and 
the prefactor of the inertial term that arises from the breaking of galilean invariance in these 
models. We show that this prefactor can be set to unity in the generalized model, therby 
effectively restoring galilean invariance. Moreover, we derive an expression for the kinematic 
viscosity, and show that it tends to decrease with the maximum number of particles allowed 
in each direction, so that higher Reynolds numbers may be achieved. Finally, we derive 
expressions for the statistical noise and the Boltzmann entropy of these models. 

I. LATTICE GASES 

Lattice gas automata (LGA) are a class of dynamical systems in which particles move on a lattice in 
discrete time steps. If the collisions between the particles conserve mass and momentum, the coarae-jgrained 
behavior of the system can be shown to be that of a viscous fluid in the appropriate scaling limitEJ u. Used 
as an algorithm for simulating hydrodynamics, the method has the virtues of exact conservation laws, and 
of unconditional numerical stability. 

In a typical LGA, there is an association between the lattice vectors and the particles at each site. If there 
are n lattice vectors, then the state of the site is represented by n bits. Each bit represents the presence 
or absence of a particle in the corresponding direction. At each time step, a particle propagates along its 
corresponding lattice vector and then collides with other arriving particles at the new site|^. The collisions 



^Note that rest particles can be subsumed into this scheme by associating them with null lattice vectors. 
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are required to conserve particle mass and momentum. 

Relevant dimensionless quantities of a LGA are the Knudsen number, Kn, defined as the ratio of the mean- 
free path to the characteristic length scale; the Strouhal number, Sh, defined as the ratio of the mean-free 
time to the characteristic time scale; the Mach number, M, defined as the ratio of the characteristic velocity 
to the speed of sound; the Reynolds number, Jle ~ M/Kn; and the fractional variation of density from its 
average value, Ap/p- Hydrodynamic behaviom is attained in the limit as Kn and Sh go to lexo. -Viscous 
hydrodynamicsO is attained when Sh ~ Kn^ in this limit. Incompressible viscous hydrodynamicstl is then 
attained when we also have M ~ Kn so that Re ~ 0(1), and 5p/ p ^ Kn^. 

The Chapman-Enskog procedure is a perturbation expansion in the above-described asymptotic ordering. 
For a LGA whose collisions conserve mass and momentum on a lattice of sufficient symmetry (Quantified 
below), the local equilibrium distribution function can be shown to be Fermi-Dirac in natureoD. The 
Chapman-Enskogrpeocedure can then be used to compute the correction to this Fermi-Dirac distribution 
and thereby showErH that the pressure, P, and the momentum density, u, obey the following equations in 
the asymptotic limit: 

V • u = 

5^ . vu = - VP + Kp)vV 

at p 

where p is the fluid density (a constant in this limit). The analysis also yields expressions for the functions 
g{p) and ^{p), and an equation of state for P. In particular, the form of these equations, the equation 
of state, and the expression for the function g[p) depend only on the fact that mass and momentum are 
conserved - and are the only things conserved - by the collisions. The expression for the viscosity, i'{p) 
depends on the details of the collision rules used. 

Since the fluid density p is a constant in the asymptotic limit, the factors g{p) and v{p) are also constants. 
As has been noted, the latter is the fluid viscosity. The presence of the former is reflective of a breaking of 
galilean invariance, due to the fact that the lattice itself constitutes a preferred galilean frame of reference. 
For a single-phase LGA, the former factor can easily be scaled away by redefining the momentum density 
and pressure as 

U = g{p)n 

and 

V = g{p)P, 

where u and P are those measured in the simulation. For compressible flow, or for multiphase flow with 
interfaces, however, the presence of the g{p) factor is problematic, and various techniques have been proposed 
to remove it^It has been shown that this can be done by judiciously violating semi-detailed balance in the 
collision rulcQ, or by adding many rest particles at each sitelS. 

The unconditional stability of the lattice gas procedure arises from a requirement that the collisions 
satisfy a statistical reversibility condition known as semi- detailed balance (SDB). The collision process is 
fully specified by the transition matrix A{s s') which is the probability that the incoming state s will 
result in outgoing state s' . Since collisions must result in some outgoing state, conservation of probability 
requires that 

5]A(.^,s') = L (LI) 

SDB is then the condition that 

Y^Ais^s')^!. (L2) 
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(Note that the condition of detailed balance (DB), A{s s') = A{s' — > s), imphes that of SDB, but not vice 
versa; that is, SDB is a weaker condition than DB.) From SDB, it is possible to prove an iJ-theorem, from 
which foUows the unconditional stability of the lattice gas algorithm. 

An important limitation of the lattice gas procedure has to do with the statistical noise associated with 
the coarse-grained averaging that is necessary to get the hydrodynamic quantities that obey the above 
fluid equations. For n bits per site, and for coarse-grained averages over blocks of N sites, the noise is 
of order ^ 1/ \J nN . For some applications - most notably the simulation of complex fluids - a certain 
controllable amount of noise is actually desirable because it is essential to the physics; for simple fluid 
dynamics computations, on the other hand, the noise is a nuisance. 



Because of their noise and lack of galilean invariance, LGA have been replaced, by Lattice Boltzmann 
Equations (LBE) for many hydrodynamics applications of interest in recent yearsS. These methods keep 
track only of an averaged occupation number of particles in each direction at each site. Moreover, the collision 
operator most often used is a simple relaxation to a noiseless equilibrium, thereby eliminating the statistical 
fluctuations that are inherent in the LGA method. This means that in complex fluid applications fpc which 
statistical fluctuations are an essential part of the physics, they have to be reintroduced artificiallyt2l. 

For a lattice Boltzmann equation corresponding to a lattice gas with only one bit per lattice vector, this 
real-valued distribution function is bounded between zero and one. This need not be the case, however, and 
the LBE procedure allows one to tailor the equilibrium distribution-function to satisfy certain desiderata. 
Among these is the ability to demand galilean invariance [g{p) = 

At the same time, the LBE method gives up two of the principal advantages of LGA's: Due to the roundoff 
error inherent in manipulations of real numbers on a computer, it no longer maintains the conservation laws 
exactly. Moreover, LBE's are no longer unconditionally stable; indeed, they are subject to a variety of 
numerical instabilities, most of which are not well understood. 



In this paper, we investigate a simple generalization of the lattice gas concept that can be used to control 
the level of statistical fluctuations - reducing it if desired, but not necessarily eliminating it altogether - 
while maintaining the conservation laws exactly, preserving unconditional stability, and allowing for galilean 
invariance. 

The use of a single bit per each of n directions to represent the state of a given lattice site means that 
the number of particles moving along any lattice direction is either zero or one. We generalize this by 
allowing for up to L bits per direction, for a total of nL bits per site, so that the number of particles moving 
along any lattice direction can range from to 2^ — 1. The total number of states per site is then 2"^. 
Computationally, this means that the state of each direction is described by an integer of L bits; hence, the 
terminology. Integer Lattice Gas Automata (ILGA). 

To simplify the derivation of the hydrodynamic equations of an ILGA, we use the Boltzmann molecular 
chaos approximation, so that all quantities in our analysis are ensemble averaged, and we indiscriminately 
commute the application of this average with the collision process. We also assume that the particles are of 
unit mass. Denote the ensemble-averaged value of the £th bit in the ith direction by A^*'^, where < i < n—1 
and < £ < L — 1. Also, denote the lattice vector for the ith direction by c^. The distribution function for 
the total number of particles in each direction is then. 



II. LATTICE BOLTZMANN EQUATIONS 



III. INTEGER LATTICE GAS AUTOMATA 



L-l 




(III.3) 



The ensemble-averaged mass and momentum densities are then given by 
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n — 1 n—1 L — 1 



i=0 



and 



n-l 



=0 e=o 



n-l L-1 



(III.4) 



(IIL5) 



i=0 



1=0 «=0 



Let us also associate an energy Si with each particle in direction i. The ensemble-averaged energy density is 
then given by 



i-l L-l 



(IIL6) 



=0 e=o 



IV. THERMODYNAMICS OF THE INTEGER LATTICE GAS 

We first consider the thermodynamics of the integer lattice gas. The grand canonical partition function is 

^ = ^exp[-/3(£-a-P-MM)], 

{N} 

where the sum is over all possible states of the ILGA (that is, each N^{x) is summed from to 2-^ — 1), 
where (3, cx and fi are Lagrange multipliers, and where 

V n 



V n 



and 



V n 



are the total mass, momentum and energy, respectively, of all the particles on a lattice of V sites. Thus, we 
have 



V n 



^ = E 

{JV} 

V n 

= E n n (£i - a • c. - /x) iV^x)] 

{JV} X i 

= fin E « • - M) fc] = nn E (^')* 



X i k=0 



fc=0 



1 - (z^) 



where we have defined the fugacity 



= exp [-P (si ~ a-Ci - n)] 



(IV.7) 



The grand potential is then 



n = — hiz = 



^ \ 1 



1 - 



/3 



so that 



d(3 ^' ' dp 
We identify the equihbrium distribution function 

z 2^{zf 



1 - (Z»)2^ 



Fl{z) 



l-Z 1 - (z)2 



(IV.. 



which gives the mean number of particles moving in each lattice direction. Since this has a maximum of 
2^ — 1, we also define the fractional occupation number 



Fl{z) 
2^-1 



Figure]^ shows fhiz) plotted against z for several values of L. 




0.5 



1 . 5 



FIG. 1. fL{z) versus z for several values of L. The black curves represent L values from 1 to 6, with increasing 
steepness, while the gray curve is the limit as L ^ cxd. 
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In terms of the equilibrium distribution function, we have 

n + (3—^n^T—^vY,FL{z'){e,~a■c,-^l). 

i 

where T = 1//? is the temperature. It follows that 

n = {H)-oi■{v)-^i{M)-T{s), 

where we have identified the average energy 



the average momentum 



the average mass 



and the average entropy 



{m)^vJ2fUz1, 



dT 

In the expression for the entropy we have defined the function 

Sl{z) ^ \n (l -z^') + ( In (^'') - In (1 - ^) - (y^) In [z) (IV.9) 

as the entropy per lattice direction. Thus, in addition to the form for the equilibrium distribution function, 
this analysis has provided us with an expression for the entropy that is additive in the contributions from each 
lattice direction. In fact, it is straightforward to show that 5*^ ^ L In 2 in the limit of large L, corresponding 
to a dominant contribution of In 2 per bit of state. The excess 

AS'L = 5L-Lln2 (IV.IO) 

is then 0(1) in L and is plotted against the fractional occupation number fhiz) in Fig. |[ This can be 
interpreted as indicating that the bits are most random at half filling; elsewhere, the entropy is lower than 
i In 2 per bit. 
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s_L - L*ln(2) 




FIG. 2. Entropy excess ASl versus fractional occupation number /^(z). The black curves decrease with increasing 
L, and the gray curve is the limit as L ^ cxj. 



V. KINETIC-THEORETICAL TREATMENT 



As an alternative to the preceding thermodynamic treatment of the integer lattice gas, we can derive the 
principal results from a kinetic-theoretical argument. For example, to derive the form of the equilibrium 
distribution function, Eq. ([V.5), we can note that the equilibrium distribution function for each bit must 



still be Fermi-Dirac in form, since each individual bit is either occupied or not. Thus, 



N, 



1 + exp [2^P {si - a - Cj - n)] ' 



where the multipliers (3, a. and /x are determined in terms of the mass, momentum and energy densities by 
their definitions, Eqs. (III. 4), ([ITS) and (IITG). In terms of the fugacity, Eq. (IV. 7), the above may be 
written, 



1 



1 + iz'Y 



(V.ll) 
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The equilibrium distribution function for each direction is then given by Eqs. ( pL^ ) and ( [vTl| ) 

L-l 



1=0 

where we have defined the function 

L-l 



In Appendix ^ we show that this sum is equal to the closed form derived in the previous section, 

Fl{z)^j^-^^. (V.13) 
1 — z 1 — 



VI. FORM OF THE HYDRODYNAMIC EQUATIONS 

To derive the hydrodynamic equations, we first expand the equilibrium distribution function in the Mach 
number. Here and henceforth, we specialize to the case of no internal energy, so that Si — 0, and we can 
absorb the multiplier (3 into a. and Treating cc as a small quantity, the fugacity can be written 

= ^1 + a • Ci + iaa : CjCj^ = zq + z\ + z^, 

where the subscripts of 

zo = e^' 

Z\ = ZqOL ■ Ci 
^0 

Z2 = ■ CiCi 

denote the order of the Mach number expansion, and we note that zq is independent of the direction i. It 
follows that 

iV^ = Fl{z') = FL{za + z\ + zl). 

Taylor expanding, we get 

Nq = Fl{zo) + zoF'^{zo)oL ■ Ci + izo [2;oF[(zo)]' aa : CiCj. 
To proceed, we must make some assumptions about the symmetries of the lattice. We demand that 

n— 1 k 

^(g)c, -Afclfc (VI. 14) 

1=0 

for < fc < 4, where ® denotes the outer product, and where 1^ is the completely symmetric and isotropic 
tensor of rank 

In = 1 



(ll). = 
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(l2)y — % 



(l3)iifc = 



i'i-4)ijki = SijSki + SikSji + SiiSjk- 



Note that Eq. (VI. 14) defines the coefRcients for a given lattice. 
We now demand that 



i=0 



and 



u = ^ CiN^ = A2ZoF'L{za)a. 



If we now let z denote the solution to the equation 

P 



i-^FLiz), 

^0 



(VI. 15) 



it follows that the difference between z and zo is of second order in the Mach number, so that we can solve 
for /i and a. We find that 



A^zFiizY 

and that /i = In zq where zq is the solution to the equation 

-f'Ll^oj - -j 



Ao 2A0A2 [zF[{z)f 



Inserting these results into the distribution function, we find 

5 u-c, z[zF'^{z)]' 



A. 



Ao A2 2Al[zF'^{z)Y 

where, again, z is defined by Fl{z) — p/Aq. 

The inviscid part of the stress tensor is then given by 



>^ A2 z\zF'r(z)]' /A4 A2 



Ao 



Ao' 



Ao 2A2 [zF[{z 

'A2 zjzF'^jz)]' /A, A2^^2 
Ao 2A2[zF'^{z)f \A2 Ao 

= P(p,u)l2+5(p) — , 

p 



A4z[zFi{z)]' ^ 
Al [zFlizf 



(VI.16) 



where we have identified the factor that multiplies the inertial term in the Navier-Stokes equations, 

AoAizFL{z) [zF^(z)]' 



9{p) 



A\ \zF'^{z)X 



(VI.17) 
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and the equation of state, 



A2 



1 - 



(VI. 18) 



Eqs. ( |VL15| ), ( |VLT7| ) a nd (|VLT1 ) are the principal results of this secti on. Eq . ( |VI.15|) giv es p in terms 
of the parameter z. Eq. (VI. 17) then gives g in terms of z, so that Eqs. (| VI.15[) an d (VI. 17) are a pair of 
parametric algebraic equations for g in terms of the density p. Finally, Eq. ( VI. 18 ) gives the equation of 
state for P in terms of p and u. The c oefficie nts Aj that appear in these equations are given in terms of the 
lattice vectors by the conditions, Eq. (VI. 14). 



VII. EXAMPLE: BRAVAIS LATTICE 



As a concrete example of this formalism, we consider the ca se of a regular Bravais lattice. Examples of such, 
lattices with the requisite symmetry conditions, Eq. ( VI.l^ , are the triangular lattice in two dimensionsEJ 
and the face-centered hypercubic lattice in four dimensionsH. In addition to the n directions corresponding 
to unit-speed particles, we include null lattice vectors to accomodate rest particles. In this situation. 



^0 



n + fir 



A2 = — 

D 



and 



A. 



D{D + iy 



where D is the number of dimensions. Inserting these into Eqs. ( VI. 15 ) through ( Vl.lg ), we find 

p = in + Ur) Fl{z), 



9{P) 



D 



D + 2 



Gl(z), 



and 



P{p,u) 



Dnr\ u 

-TT- gyp)— 

2n J p 



21 



Here we have defined the function, 



z,fL{z) [zfl{z)\ 

which we plot against the fractional occupation number, 

P 



fUz) 



^ Fl{z) 



(VII. 19) 



for several different values of L in Fig. |^. For L = 1 it is a straightforward exercise to show that we recover 
the well known resuliij. 



Gi(z) = 



1-2/ 
1-/' 
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which decreases monotonicany from unity at / = 0, to zero at half-fiUing (/ = 1/2), after which it becomes 
negative. For L > 1, we see that this decrease is no longer monotonic, since the slope at the origin, g'{0), 
is positive. Thus, for L > 1, the function g has a maximum for some < / < 1/2. The location of 
this maximum approaches / = as L — > oo. (The limit of infinite integers, i.e., L — > oo is discussed in 
Appendix ^ and is shown as a shaded curve in Fig. ^.) 



G_L 




0.1 0.2 0.3 0.4 0.5 ~" 



FIG. 3. Gl versus / for several values of L. The black curves represent L values from 1 to 6, increasing upward, 
while the gray curve is the limit as L ^ oo. 

Galilean invariance is achieved when g = 1, or 

If the quantity (1 + 2/ D)n/(n + n^) is greater than the maximum value of Gl, then galilean invariance is 
impossible for those values of D, n, rir, and L; if it is less than this maximum, then there are two densities 
at which galilean invariance is achieved. Some of these values are tabulated for the FHP and FCHC lattice 
gases in Fig. ^ 
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FHP Lattice Gas (D = 2, n ^ 6) 
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Low-Density Root High-Density Root 
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0.212636 


4 
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0.276535 




FCHC Lattice Gas (£> = 4, n = 24) 
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FIG. 4. Values of / G (0, 1/2) such that g = l 



VIII. VISCOSITY 

To compute the viscosity of a ILGA in the Boltzmann molecular chaos approximation^, we consider its 
ensemble-averaged collision operator, fi*'^. This quantity is the ensemble average of the increase in bit £ in 
direction i due to collisions. It is given by 

s.s' 

where V{s) is the probability that the incoming state is s, A{s — > s') is the probability that the collision 
process takes incoming state s to outgoing state s', and s*'^ is the value of bit £ in direction i in incoming 
state s (and likewise for outgoing state s'). In the Boltzmann approximation, the probability of a state s is 
the product of the corresponding fractional occupation numbers, or their complements, 

n— 1 L — 1 ^_ k\j' 

^(^) ^ n n (^'''0 {i-N''^=')~ . 

k'=Of=0 
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To get the total increase of particles in direction i, we take the sum 

L-l n-1 L-1 k',j' 

EE 2'^^''' = Y,A{s~^ s'){s'' - s') Y[ Y[ (A^'-'-^y (l - iV'=''- 

£=0 s,s' k'=Of=0 



l-s" 



where 



L-l 
1=0 



is the total number of particles in direction i in state s (and likewise for s'). 

To compute the viscosity, we must form the Jacobian matrix of the collision operator. Direct calculation 
yields 



We would like to evaluate this Jacobian at the equilibrium given by Eq. (V.ll), 



N, 



l+(z'=)-2^' 

Taking the derivative of this equation with respect to the fugacity, 



dN, 



k\-V-\ 2^A^n ( 1 - iV, 



we can use the chain rule to get the integer version of the Jacobian of the collision operator at equilibrium, 





•J i 



L-l 



dNf'/dz^ 



k _ jsfk\ 



(VIII.20) 



In fact, we need this result only in the limit of zero Mach number, so we can use the lowest order expression 
for the fugacity, = z (see Sec. VI), which is independent of the index k. We find that the zero Mach 
number limit of the Boltzmann probability of state s is given by 



n-l L-l ^fc',/ 

^o(-) = n n K'O (i 

fc'=0 j'=0 



N, 



Hi 



+ z- 



'2i' 



pAs) 



l + z- 



-2i' 



where 



k,3 



k=Q 



is the total number of populated bits in the jth binary digit. It follows that 
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L-l 



where 



L-l 

p{s) ^ 



1 - z 

1 -Z2^ 



,P(s) 



is the total number of particles present in state s. 



Inserting this result into the expression, Eq. (VIII. 20), for the collision operator, we obtain 



zF'^{z) \\-z 



1- z 



J2 Ms s') {s" - s') [s'^ - Fl{z)] zP^'K 



s,s' 



As a consequence of conservation of probability, Eq. (p]), and semidetailed balance, Eq. (p^), it follows that 
the second term in square brackets vanishes, so we finally get 



1 



zF^^iz) yi-z^ 



1 - z 



At first order in Knudsen number, the kinetic equation iM 

c, • VN^ = J)N(, 

where there is an understood summation over j. The only part of th e left-h and side that contributes to the 
viscosity comes from the second term on the right-hand side of Eq. ( VI. 16 ), whence 

Now, J is a singular matrix; it has a null eigenvector corresponding to each hydrodynamic mode of the 
system. These null eigenvectors span what we shall call the hydrodynamic subspace of the system. The 
complement of this subspace is called the kinetic subspace, and is spanned by the kinetic modes with nonzero 
(negative) eigenvalue. If we restrict our attention to the kinetic subspace, then we can form the pseudoinverse 
of J, denoted by J~^, in terms of which we may write 



A 



'j ,CjCj : Vu. 



The conservation law for momentum then contains the term 
(c»c, • VNl + ic,c,c, : VViV^ + • • •) 



^ ^ C-^C^ ^t/ ^ jCjCj 2 ^ ^ C^C^C^C^ 



(Vu) 



We note that J ^ is diagonalized and degenerate in the subspace spanned by the n outer products of the 
lattice vectors with themselves; that is 



(VIII.21) 



where A is a scalar, whence the above term in the momentum conservation equation becomes 
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V • 



-l(-A + i)l4:(Vu) 
A2 



= V 



4^(-A + i)Vu 
^2 



from which wc identify the kinematic viscosity, 



The quantity A is then determined by taking the double spatial dot product of CiCi on both sides of 
Eq. (VIII. 21), and summing over i to get 



whence 



-1 / 1-z 



= -A5:J^(c 



where / — fhiz) determines the parapaeter z in terms of the fractional occupation number. This result is 
easily seen to reduce to that of H6noi£3 when L = 1. 

We computed the viscosity of an L = 2 lattice gas in two dimensions {D = 2) by measuring the decay 
of a sheas, wave in periodic geometry. We used a lattice of size 512 x 512 on a CAM-8 Cellular- Automata 
MachineEJ. The probabilistic collision procedure used obeyed semi-detailed balance, with each outgoing state 
allowed by the conservation laws sampled with equal probability. Fig. || shows the decay of the shear wave 
amplitude to be exponential in nature, as is appropriate for Navier-Stokes evolution. The time constant of 
the exponential then determines the viscosity, which is plotted as a function of density in Fig. along with 
the curve predicted by the theory given above. 
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FIG. 5. Time decay of shear- wave amplitude 
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Fractional Occupation 

FIG. 6. Viscosity versus / for L = 2. 

While the agreement with theory is good at intermediate values of the fractional occupation number near 
half filling, we note that it is seriously in error at low (and high) fractional occupation numbers. At present, 
we attribute this discrepency to deviations fram the Boltzmann molecular chaos approximation, and we plan 
to investigate them using kinetic ring theoryo in a forthcoming publication. 

IX. STATISTICAL NOISE 

Finally, we consider the statistical noise of the ILGA model. With the maximum number of particles per 
direction increasing as 2^, one might naively expect the noise level to decrease with L as 1/V2^ ~ 2~^/^. 
Unfortunately, as we shall show, this expectation is not realized, due to the extremely narrow dynamic range 
of the fugacity for large L. This is best seen in Fig. in which the effective width of the function fL{z) near 
z = i decreases like 2~^, making for a subtle limiting process that is discussed in Appendix ^ 

Let n*'^(x, t) be the precise value of bit £ in direction i at lattice site x at time t. The ensemble-average 
of this quantity is N"^'^, as used in the text of the paper. The mean number of particles in a (space-time) 
block of N sites is then 

N n L-1 

(x,t) i e=o 

where the angle brackets denote the ensemble average. 

The mean square of the number of particles in this block of sites is then 

N N n n L-1 L-1 

^^-E E EEEE2^"'(«^''(-'*K''(-''0). 

(x,t) (x',t') i i' 1=0 i'=0 

The bits are either zero or one, and in the Boltzmann molecular chaos approximation different bits are 
uncorrelated. It follows that 
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(n^'^(x,t)n^''^'(x',i')) = (n^''(x,t)) (n''-^' (x', t')) + <5x,x'<5m'5»,.'<5a^' (n^''(x,t)) (l - (n^'^(x,t))) 
whence 



^2 = + niV ^ 



= Tl+nNzFl{z). 



It follows that the standard deviation of the number of particles in the block is \/ T^ — T\ ■ To define a 
fractional noise, we could divide this by the mean number of particles, J-\, but it preserves particle-hole 
symmetry if we instead divide it by the square root of the product of the mean number of particles and the 
mean number of holes, thus 



1 



^1 {nN (2^ ~ 1) - Ti\ ^n7V(2i-l) V /i(^) [1 " h{z)\ ' 



(IX.22) 



This appears to decrease exponentially with L, but it must be noted that the logarithmic derivative of /l(2) 
goes as 2^ at z = i. Since, for fixed fractional occupation number f^, z tends to ^ as i tends to infinity, 
we see that AjF is order unity in L. Thus, the fractional noise does decrease with L, but not as rapidly as 
one might hope. It is plotted for several different values of L in Fig. |^. 
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Fractional Occupation Number 

FIG. 7. ^Tl versus / for several values of L. The black curves represent L values from 1 to 6, increasing downward, 
while the gray curve is the limit as L ^ oo. 
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X. SAMPLING PROCEDURE 



Finally, we consider some practical considerations concerning the computer implementation of the ILGA 
model. Since each site has nL bits, and therefore 2"^ possible states, and since the most popular lattices 
with the requisite isotropy properties have n = 6 and n = 24, it is clear that the brute-force approach in 
which a lookup table is used to store the collision outcome states will not be feasible for L much greater 
than unity. 

For this reason, we propose another sampling scheme for the outgoing states. Though the method we 
propose is completely general, we illustrate it for the two dimensional integer lattice gas on a triangular grid 
(n = 6). Let n be an integer- valued column n- vector whose components are the particle occupation numbers 
in each of the six directions. 

Let us suppose that the mass and the two components of momentum are the only conserved quantities. 
Since these conserved quantitites are linear in the particle occupation numbers, each of them correspond to 
a row vector, whose inner product with n yields the conserved quantity in question. Thus, corresponding to 
the mass we have the row vector 

qi = ( 1 1 1 1 1 1 ) , 
corresponding to the a;-momentum (multiplied by a factor of 2), we have 

q2 = ( 2 1 -1 -2 -1 1 ) , 

and corresponding to the y- momentum (multiplied by a factor of 2/\/3), we have 

qa = (0 1 1 -1 -1). 

In fact, these row vectors are precisely the hydrodynamic eigenvectors, mentioned in our derivation of the 
viscosity; that is, 

A(qi)' = J'M' = J'M' = 0. 

It is clear that these can always be chosen to be mutually orthogonal, without loss of generality. Using the 
Gram-Schmidt procedure, it is then possible to find three vectors spanning the kinetic subspace, orthogonal 
to the above; e.g., 

q4 = (2 -1 -1 2 -1 -1), 
q5 = (1 -1 1 -1 1 -1), 

and 

qe = ( 1 -1 1 -1 ) . 

Now the collision process takes state n to state n'. Since it cannot change the values of the conserved 
quantities, it follows that the difference between n' and n must be a linear combination of kinetic eigenvectors. 
That is, 

n' = n + 0:4 ql + a^ql^ + aeqj , 

where the a's are integer constants, and where the superscript T denotes "transpose." Thus, writing out 
components, we have 

/ ni + 2a4 + a5 \ 

n2 — a4 — + ag 
^, _ ns — a4 + as — 

Ui + 2a4 — Q:5 

n5 — a4 + as + ae 
\ ne - Q!4 - Q!5 - Q!6 / 
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Since the components of n' must all be between and 2^ — 1, inclusive, we derive the following six inequality 
constraints: 

< ni + 2a4 + as < 2^ - 1 
< n2 — — + ae < 2^ — 1 
< n3 — + " ae < 2^ ~ 1 
< ^4 + 2a4 - as < 2^ - 1 
< ns - a4 + as + ae < 2^ - 1 
< ?^6 ~ a4 — as — ag < 2^ — 1. 

These inequality constraints define a polytope in the three dimensional space of allowed values of a4, as and 
ag. We know that this polytope must exist and contain the origin in that space, since a4 = as = ae = 0, 
corresponding to the "trivial collision" in which the occupation numbers do not change their values, will 
always satisfy the constraints. 

The collision process is then specified by a strategy for sampling points from this polytope. One viable 
strategy which certainly satisfies semidetailed balance is to sample the points within this polytope uniformly. 
It is possible, though tedious, to derive a closed-form algorithm to do this, based on the above constraints. 
Alternatively, with some loss of efficiency, one can simply bound the polytope sffld use a rejection sampling 
scheme. Details of this procedure will be provided in a forthcoming publicationllil. 



XI. CONCLUSIONS 



We have generalized the hydrodynamic lattice gas model to include integer numbers of particles moving 
in each direction at each site. We have presented the thermodynamics and kinetic theory of this generalized 
Integer Lattice Gas (ILGA) model, including closed-form (or parametric algebraic) equations for the equi- 
librium distribution function, the entropy, the equation of state, the non-galilean factor in the inertial term 
of the fluid equations, and the statistical noise. We have thereby shown that the ILGA model allows for the 
attainment of galilean invariance, and a reduction in the kinematic viscosity and the statistical noise. In 
future publications, we shall show that this generalization also allows for more straightforward inclusion of 
interparticle interactions than the usual binary model. 
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APPENDIX A: CLOSED-FORM EXPRESSION FOR Fl{Z) 



In this Appendix, we prove Eq. ( V.13 ), where Fl{z) is defined by Eq. ( V.12| ). Using mathematical 
induction, we first note that the statement of the theorem is true for L ~ \: 



z z[l — z) 



z' 



Z 



z{l + z) 2z2 



z 



1 + Z 1 - z2 

2z2 



1 — z^ 1 — 1 — z 1 — 

Next, we assume the truth of the statement ioi L = K: 

K-l „f 



Fk{z) 



z 



2«Z 



+ z 



-2* 



1 - Z 1 - z2 



It follows that 



K 



2^ 2^ 

Fk+i{z) = J2 T-—¥ = + 1 

^1 + z 1 



z 



z 



2^z2' 



1 - Z 1 - z2 



l + z- 



1 - z 



1 - z- 



TT + 



1 + z- 



2K+1 



1 - Z 



1- (z-2'^)^ 



1 - z 1 - z- 



_2A- + l 



1 — Z 1 — Z^ 

and we have thereby proven the theorem for all K. 

Alternatively, we may simply note that the summation can be written in the telescoping form 



1=0 l=Q \ 



L-1 



1-z 



2«+l I ' 



from which the result follows immediately. 
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APPENDIX B: THE INFINITE INTEGER LIMIT 



To consider the limit of infinite integers, L — > oo, we first note that the fractional occupation number, 

Fl{z) 1 / Z 2^22"" \ 



2^-1 2^ - 1 I 1 - 1 _ ^2^- I ' 



(B.l) 



has the limiting behavior 



lim/L(z) = 



lim fL{z) = \ 



lim /l(z) = 1 

Z — 'OO 

for all L] here we have used L'Hopital's rule to establish the result for z ^ 1. Referring to Figure |l|, we note 
that the function fhiz) becomes increasingly like a step at z = 1 as L — > cx). To verify this, we note that 
the width of the gradient there can be estimated by 



lim 



fdz) 



if'Jz) 2^ + 1' 

which clearly goes to zero as L — > oo; once again we have used L'Hopital's rule to establish this result. 

The approach to a step function means that the entire range of fractional occupation numbers is 
parametrized by values of z within order 2^^ from 1, as L — > oo. That being the case, we write 



1 



y_ 

2^' 



(B.2) 



where y is a new par amet er of order unity. Note that the fractional occupation number is exactly 1/2 when 
y = 0. Inserting Eq. (B.2) into Eq. (B.l), we can now take the limit as i — > oo to get 



- e-y 



(B.3) 



Next, inserting Eq. ( |B.2| ) into Eqs. ([V.9) and ( IV.10| ), and taking the limit as L ^ oo, we find the entropy 
excess. 



(B.4) 



Eqs. (B.3) and (B.4) constitute parametric algebraic equations, with parameter y, yielding AiSl as a function 
of the fractional occupation number as L ^ oo. These equations were used to produce the shaded curve 
in Fig. |. 



Likewise, inserting Eq. (B.2) into Eqs. ([X.22), and taking the limit as L — > oo, we find the fractional 



lim + ^ 



J/2 — 2 cosh y + 2 



y"^ — 2y sinh y + 2 cosh y — 2 



(B.5) 



Eqs. ( p3.3D and (B.5) constitute parametric algebraic equations, with parameter y, yielding as a function 
of the fractional occupation number as L — > oo. These equations were used to produce the shaded curve 
in Fig. 0. 
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Finally, inserting Eq. (B.2) into Eqs. (VII. 1£), and taking the limit as L — > oo, we find the Gl factor for 
a Bravais lattice, 



lim Gz, 



2 + + 2y - s) e" + i^y^ - Gy + I2) e^" + (y" _ y3 + 61; ~ s) e^" - 2 (y - 1) e'' 
[1 - + 2)e!' + e2a]2 



(B.6) 



Eqs. (B.3) and ( p.q ) constitute parametric algebraic equations, with parameter ?/, yielding Gl as a function 
of the fractional occupation number /l as L — s- cx). These equations were used to produce the shaded curve 
in Fig. I 
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